////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//---------------------------- Ebrahim Foulaadvand, 10 July 2013 -----------------                        //
//                                                                                                        //
//  The routine "AdvectionLax" solves the 1D wave equation using the Lax and Lax-Wendroff methods.        //
//  Periodic boundary condition is used. The system length L and the wave velocity c are set to one.      //
//  Initial wave shape is cosine-modulated Gaussian pulse.                                               //                         //
//                                                                                                        //
//                                                                                                        //                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>

using namespace std;

main()
{

double tau=0.004,h,L=1,c=1,sigma=0.05,lambda=0.2,k=2*M_PI/lambda,x;

int N=200,i,n,T=1000;


ofstream file0 ("initial profile n=0.plt");  //output file for the profile at timestep n=0.
ofstream file1 ("profile Lax n=20.plt");    //output file for the profile at timestep n=100.
ofstream file2 ("profile Lax n=40.plt");    //output file for the profile at timestep n=500.
ofstream file3 ("profile Lax n=100.plt");    //output file for the profile at timestep n=500.
ofstream file4 ("profile LW n=20.plt");     //output file for the profile at timestep n=500.
ofstream file5 ("profile LW n=40.plt");    //output file for the profile at timestep n=500.
ofstream file6 ("profile LW n=100.plt");    //output file for the profile at timestep n=500.



vector<double> u(N+1,0),unew(N+1,0),uL(N+1,0),uLnew(N+1,0),uLW(N+1,0),uLWnew(N+1,0);

 // Arrays "u" and "unew" store the current and updated values
 // of the scaler function u(x,t) at grid points. "L" stands for Lax and "LW" stands for Lax-Wendroff.

 h=L/N; cout<<"h= "<<h<<endl; cout<<"tau= "<<tau<<endl; //getch();


for(i=0;i<=N;i++){

x=-0.5*L+i*h; u[i]=cos(k*x)*exp(-0.5*(x/sigma)*(x/sigma));

uL[i]=cos(k*x)*exp(-0.5*(x/sigma)*(x/sigma)); uLW[i]=cos(k*x)*exp(-0.5*(x/sigma)*(x/sigma));

file0<<(-0.5*L+i*h)<<"  "<<u[i]<<endl; }

//storing the initial condition data in the data file.



for (n=1;n<=T;n++){ //cout<<"n= "<<n<<endl; //getch();


// -------------------       FTCS method -------------------------------------

for(i=1;i<N;i++){ unew[i]=u[i]-(c*tau/(2*h))*( u[i+1]-u[i-1] ); }

unew[0]=u[0]-(c*tau/(2*h))*( u[1]-u[N] ); // periodic condition has been used: u[-1] is replaced with u[N]

unew[N]=u[N]-(c*tau/(2*h))*( u[0]-u[N-1] ); // periodic condition has been used: u[N+1] is replaced with u[0]


//if (n==10)  for(i=0;i<=N;i++){ file1<<(-0.5*L+i*h)<<"  "<<unew[i]<<endl; } //storing in the data file.

//if (n==20) for(i=0;i<=N;i++){ file2<<(-0.5*L+i*h)<<"  "<<unew[i]<<endl; } //storing in the data file.


//-----------------------------     Lax Method      ------------------------------------

for(i=1;i<N;i++){ uLnew[i]=0.5*(uL[i+1]+uL[i-1])-(c*tau/(2*h))*( uL[i+1]-uL[i-1]); }

uLnew[0]=0.5*(uL[1]+uL[N])-(c*tau/(2*h))*( uL[1]-uL[N]);// periodic condition has been used: uL[-1] is replaced with uL[N]

uLnew[N]=0.5*(uL[0]+uL[N-1])-(c*tau/(2*h))*( uL[0]-uL[N-1]);// periodic condition has been used: uL[N+1] is replaced with uL[0]


if (n==20)  for(i=0;i<=N;i++){ file1<<(-0.5*L+i*h)<<"  "<<uLnew[i]<<endl; } //storing in the data file.

if (n==40)  for(i=0;i<=N;i++){ file2<<(-0.5*L+i*h)<<"  "<<uLnew[i]<<endl; } //storing in the data file.

if (n==100)  for(i=0;i<=N;i++){ file3<<(-0.5*L+i*h)<<"  "<<uLnew[i]<<endl; } //storing in the data file.



//-----------------------------     Lax-Wendroff Method      ------------------------------------

for(i=1;i<N;i++){ uLWnew[i]=uLW[i]-(c*tau/(2*h))*(uLW[i+1]-uLW[i-1])+(c*c*tau*tau/(2*h*h))*(uLW[i+1]+uLW[i-1]-2*uLW[i]); }

uLWnew[0]=uLW[0]-(c*tau/(2*h))*(uLW[1]-uLW[N])+(c*c*tau*tau/(2*h*h))*( uLW[1]+uLW[N]-2*uLW[0]);

uLWnew[N]=uLW[N]-(c*tau/(2*h))*(uLW[0]-uLW[N-1])+(c*c*tau*tau/(2*h*h))*( uLW[0]+uLW[N-1]-2*uLW[N]);

if (n==20) for(i=0;i<=N;i++){ file4<<(-0.5*L+i*h)<<"  "<<uLWnew[i]<<endl; } //storing in the data file.

if (n==40) for(i=0;i<=N;i++){ file5<<(-0.5*L+i*h)<<"  "<<uLWnew[i]<<endl; } //storing in the data file.

if (n==100) for(i=0;i<=N;i++){ file6<<(-0.5*L+i*h)<<"  "<<uLWnew[i]<<endl; } //storing in the data file.



for(i=0;i<=N;i++){ u[i]=unew[i]; uL[i]=uLnew[i];  uLW[i]=uLWnew[i];  }


}


cout<<"end." <<endl; getch();

}